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ABSTRACT 


Sampled-data  control  systems  are  coming  of  age.  Many  methods 
have  been  formulated  by  which  both  continuous  and  sampled-data 
control  systems  may  be  optimized.  Here  one  method  of  optimum 
control  design  is  presented  using  the  approach  of  dynamic  pro¬ 
gramming,  A  program  to  generate  the  state  transition  matrices 
from  the  system  dynamics  was  developed,  A  program  for  minimizing 
rather  general  cost  functions  was  written  and  examples  are  pre¬ 
sented  utilizing  the  results • 
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CHAPTER  I 


INTRODUCTION 


In  general,  any  linear  control  system  may  be  expressed  by  the  fol¬ 
lowing  state  equation: 


=  Fx  +  Du 


(1) 


where  represents  the  system  variables  and  u^  are  the  control  vari¬ 
ables,  F  and  D  are  nxn  and  nxm  matrices  defining  the  dynamical  rela¬ 
tionship  of  the  states  and  controls.  We  will  sample  the  states  at  dis¬ 
crete  intervals  and  at  a  sampling  rate  with  period  T,  A  recursive 
relationship  from  the  solution  of  the  above  differential  equation  may 
be  established: 

x_(k+l)  =  e^x_(k)  +  x^Du_(k)di  (2) 


where 


Jn2  rp  2 

eFT=  I  +  FT  +  - -  +* 

2 : 


•+ 


pH  ^n 


(3) 


Equation  (2)  may  be  written: 

x(k+l)  =  $jc(k)  +  Au(k)  (4) 


where 


(5) 


and 


F(T-t) 
e  Ddr 


(6) 


As  an  example  of  the  above  equations,  the  open-loop  system  of 
Figure  1-1  will  be  considered. 
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Figure  1-1;  Block  diagram  of  an  open  loop  control  system. 

By  inspection,  the  following  state  equation  defines  the  example  sys¬ 


tem: 


0 

1 

i _ 

IV 

X  + 

p 

o 

_1_ 

(7) 


The  <J>  matrix  as  defined  by  (5)  has  only  two  non-zero  terms  in  the 
Taylor  expansion  for  the  exponential  function: 


♦  = 


[3 


*f  3- 


(8) 


or 


1  T 


4  = 


(9) 


(0  1J 

In  a  like  manner,  the  convolution  integral  of  (6)  may  be  solved 
for  A  : 


A  = 


r 

r 

1  T-t 

0 

l 

jP  i  _ 

1_ 

dT 


(10) 


giving 


A  = 


T2 


% 


(ID 


The  system  difference  equation  may  now  be  written  as  a  function 
of  the  sampling  period  T: 


x(k+l)  = 


1  T 

0  1 


x(k)  + 


u(k) 


(12) 


The  calculation  of  the  transition  matrix,  and  the  external 
forcing  function’s  impulse  response  matrix,  A,  from  the  matrices  F  and  D, 
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of  the  process  dynamics  was  accomplished  with  the  digital  computer.  This 
program  required  an  input  of  F  and  the  sampling  rate  T,  The  amount  of  ac¬ 
curacy  desired  may  he  varied  by  the  user  by  inputxng  the  least  significant 
figure  desired.  This  test  point  for  accuracy  results  from  the  fact  that 
the  sum  of  all  terms  in  a  Taylor  series  following  the  test  term  is  less 
in  magnitude  than  the  test  term.  Figure  1-2  is  the  flow  chart  of  the 
program.  The  program  has  been  incorporated  into  other  programs  written 
and  presented  in  this  paper. 
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STOP 


Figure  1-2:  Flow  Chart  for  PROGRAM  PHIDEL 
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The  optimal  control  of  a  system  depends  upon  a  valid  determination 
of  an  optimal  control  policy  which  will  maximize  or  minimize  a  perform¬ 
ance  criteria.  This  control  policy  may  be  generated  from  criteria  in¬ 
volving  time*  system  state  variables,  noise,  control  effort,  etc.  In 
some  cases,  this  control  policy  may  be  generated  by  active  networks  or 
filters  incorporated  within  the  system,  or  a  digital  computer  may  be 
used  to  generate  the  control.  To  accomplish  this,  all  dynamics  of  the 
system  must  be  fed  back  to  the  digital  computer  as  indicated  in  Figure  1-3. 


PLANT 


Figure  1-3:  Multivariable  Digital  Control  System 
A  general  approach  to  the  optimization  of  any  dynamic  system  was 
developed  through  the  efforts  of  R.  E.  Bellman  and  has  been  called 
,Tdynamic  programming" .  As  an  example  of  this  method,  in  a  missile  prob¬ 
lem,  one  might  choose  a  proper  cost  function  (terminal  CEP)  and  minimize 
this  cost  throughout  the  trajectory.  The  optimum  trajectory,  according 
to  the  chosen  cost  function,  would  be  computed. 

Dynamic  programming  appears  to  be  an  excellent  method  for  determ¬ 
ining  an  optimal  control  process  because  of  its  adaptability  to  use  on 
a  digital  computer.  Throughout  the  remainder  of  this  paper,  this  tech¬ 
nique  will  be  utilized  in  formulating  general  programs  to  synthesize 
optimum  controllers. 
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CHAPTER 


II 


PROGRAM  OPCON1 


Again  let  us  cons Ider  the  system: 

x  =  F x  +  D  u  (1 f ) 

whose  discrete  solution  is 


x(k+l)  =  <f>x(k)  +  Ap(k) 


C4* } 


where  we  now  take  the  case  where  u(k)  Is  a  scalar  control  and  is 
a  column  vector  of  system  impulse  responses  due  to  control  inputs. 

The  recursion  formulas  will  be  derived  for  this  system  minimizing  the 
sum  of  the  quadratic  terms  as  described  by  the  matrix  Q  over  an  N  stage 
process.  The  Q  matrix  is  specified  by  the  requirements  of  the  user  and 
the  application  of  solid  engineering  judgment.  This  cost  function  may 
be  written: 


N  T 

J(N)=  minimum  l  xi(k)Qx_(k) 
k=l 


(13) 


By  starting  with  the  last  stage  first,  equation  (7)  becomes: 

(  rp  M”  1  rT  1 

J(N)=  minimum  ^_xA  (N)Qx(N)+  ^ l  x  (k)Qx(k)j 


or 

J(N)=  minimum  £jt^(N)Qx(N)^  +  J(N-l)  (14) 

It  will  now  be  noted  that  only  the  _xT(N)Qx(N)  term  contains  the  control 
policy  u(N-l) ;  therefore,  to  minimize  this  function  over  the  last  stage, 
we  need  only  to  find: 
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3(xl(N)Qx(N)) 
3  u(N-l) 


=  0 


(15) 


By  substitution  of  equation  (4)  into  xT(N)Qx(N): 

4 

xT(N)Qx(N)=  [^x(N-1)+4u(N-1)]  TQ  [j>x(N-l)+Au(N-l)] 


and 


or 


?uHW>  =  2^T<^  u(N-1)=0 


T 

A  Q<t> 


^  2L(N-1) 


(16) 


(17) 


(18) 


Equation  (18)  now  represents  the  optimum  control  policy  u(N-l),  It 
will  be  noted  that  we  are  dividing  by  a  multiplication  of  three  matrices 
this  'is  a  valid  operation  because  a  row  matrix  times  a  square  matrix 
times  a  column  matrix  will  always  yield  a  scalar  number  if  the  multi¬ 
plication  is  valid.  The  multiplication  is  valid  if  the  dimensions  of 
the  respective  matrices  are  proper;  and  in  this  case,  it  will  always  be 
valid. 

By  proceeding  back  one  more  stage,  the  control  policy  u(N~2)  may 
be  optimized, 

J(N}=  minimum  (N)Qx(N)  +x^  (N-l)Qx(N-l)^  +J(N-2)  (19) 

Again  by  direct  substitution,  equation  (19)  becomes: 


J(N)-  minimum  \  ^  }{<f>£(N-2)+Au(N-2)  nTq 

(  L  a  Qa  J 


I{^Af^}{fe(N-2)+Au(N-2)}] 

L  A  J.  ( 

+  |$x(N-2)+£u(N-2)J  Tq  ^x(N-2)+Au(N-2)J  ^  +J  (N-2) 
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Minimizing  with  respect  to  u(N-2) ,  remembering  that  J(N-2)  is  not  a 


function  of  u(N-2);  and 


0: 


“  2ATQ<J,x(N-2)+2ATQ4u(N-2) 

+24TH  — ‘--f  4>TxT(N-2)Q{<fi  -  ^  ^  }A 

ArQA  A^A 

,  6ATQ<f>  T  Mty 

+2  A^OA  }- Q{<|)"  ATQA  }A  u(N-2)=0  (21) 


Let 


h_  AAT(Q+pn^ 
L  AT(Q+Po)A 


therefore 


<f> ,  where  P 

o 


0 


citJ(N)  „  „ 

^u(N-2)  =  A  Q<J>ac CN— 2)+A±^iQ  tj^Ax  (H-2) 


(22) 


+ATQAu(N-2)+AT^Qij/]LAu(N-2)=  0 

(23) 

or 

..  u(N  2)=  — Jr  l~~ — -  x (N-2 ) 

A  (Q+^Q^a  " 

(24) 

Let 

P  =  ,Tq* 

1  1  1 

(25) 

T 

,  ,  A  (Q+P-. )  <p 

u(N-2)=  -  "  1  x (N-2) 

AT(Q+Pi  ) A 

(26) 

The  optimum  control  policy  for  the  (N-2)  stage  has  been  found  represent¬ 
ed  by  equation  (26)* 

To  further  the  establishment  of  recursion  formulas  *  the  next  stage 
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back  must  be  considered  and  minimized  with  respect  to  u(N-3)*  The  fol¬ 
lowing  symbology  will  be  utilized  in  this  stage  of  the  development: 


and 


A  =  |^x(N-3)+Au(N-3)] 

A  r  f  A  C  Q+P  -i )  — i 

L1  1  \  at(q+p1).a'J 


The  cost  function  then  becomes: 


(27) 


(28) 


J(N)  =  minimum  {ATQA+BTQB}+J (N-3) 

(29) 

Let 

♦,  -  [i-  Mt(Q+Pi)  1* 

At(Q+Pi)aJ 

(30) 

Therefore 

J (N)=  min imum { ATQA+  TQ  ^+J(N"3) 

(31) 

Minimizing 

the  cost  function  with  respect  to  u(N-3): 

3u(N-3)  =  2ATQte(N-3)+2ATQAu(N-3) 

+2ATi|>„T(|1q^  ^  (jix(N-3) 

"2112 

+2ATV'2T,*'1TQ^1^2— u  ^N-3>  =  0 

(32) 

Solving  for 

u(N-3): 

„rN  -»v_  A^HAV ’f' 

u(N-3)-  -  21  1  ^  xfN--n 

^QA+fiT^  ^  ^  A 

(33) 
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giving 


Let 


u(N-3)  -  -  - 


at(q+*2^q^2h 

£t(Q*^2*iQ*iV£ 


x(N-3) 


p  =  ifjTlNTQaj  iL 

2  v2vl^lv2 


therefore 


u(N-3) 


at(q+p2)  * 
hT(Q+P2)  a 


x(N-3) 


(34) 


(35) 


(36) 


It  is  believed  that  an  iterative  relationship  may  now  be  shown  without 
regressing  to  another  stage.  The  following  equations  may  readily  be 
shown  to  be  true  by  an  inductive  proof. 


AAT(Q+pk)  ' 

AT(Q+Pk)  A; 


(37) 


where 


or 


P  = 


♦J*Tk-r",2*iQ*i*2"-Vi*k 


(38) 


pk  =  KP  *  ,  P  =  o 

K  k  k-1  k  ° 


and 


u(N-k) 


=  _  A  (Q+Pk.i) 


AT(Q+Pk_1)  A 


x  (N-k) 


(39) 


(40) 


x(k+l)  =  (j>x_(k)+Au(k) 


(41) 


To  illustrate  the  use  of  the  above  equations,  the  system  presented 
in  Chapter  I  will  be  considered  with  the  sampling  period  equal  to  one 


second: 

i 

1 

V 

r 

1 

o' 

♦  = 

A  = 

Q  = 

0 

1 

.1. 

0 

0_ 
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To  find  u(k)j  equations  37,  39,  and  40  must  be  solved: 


-  MTQ4>  = 
atqa 

and 

?i  =  »  0 

and 

,  .  aTqo, 

u(k)=  -  == —  x(k)  =  -2x-(k)-2x,(k) 
A  Q  A 


(42) 


(43) 


(44) 


By  substitution  into  equation  (44)  and  solving  for  u(0),  we  find 
^ 1 ( 1 ) =0  and  X2(l)=  “2*  The  following  phase  plane  (Figure  II-l)  results : 
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1  . 
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/' 

V 

Figure  II-l:  Phase  Plane  For  Example  Problem 

This  clearly  indicates  that  a  limit  cycle  results  and  that  the  cost  func 
tion,  as  determined  with  the  use  of  this  weighting  matrix,  yields  a  rela 
tive  minimum  cost  of  one  unit- 

Using  the  derived  recursion  formulas  37,  39,  40  and  41,  PROGRAM 
0PC0N1  was  written-  Figures  II- 2  through  11-10  represent  abbreviated 
flow  charts  of  the  program  and  the  subroutines  used  in  the  program-  It 
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must  be  remembered  that  the  arithmetic  operations  represented  in  the 
flow  charts  are  all  matrix  operations.  As  can  be  seen  from  Figure  II-2, 
the  required  inputs  to  the  program  must  be  read  from  data  cards.  The 
following  table  (Table  II-l)  indicates  the  order  that  the  data  cards  are 
read*  the  format  required  on  these  cards*  and  the  definition  of  the  vari¬ 
able  concerned.  The  primary  output  of  0FC0N1  is  the  value  of  each  state 
variable,  the  control  signal,  and  the  cost  of  each  stage.  A  complete  pro¬ 
gram  listing  may  be  found  in  the  appendix. 


Variable 

Definition 

Format  of 

Data  Card 

M 

Number  of  stages  plus 
one 

15 

N 

Order  of  system 

15 

A 

Column  vector  of 
system  impulse 
responses  due  to 
control  inputs 

8F10.6 

(read  from  top  of 
column  to  bottom) 

State  transition 
matrix 

8F10 , 6 

(read  in  order  by 
columns) 

Q 

Weighting  matrix  of 
the  cost  function 

8F10 , 6 

(read  in  order  by 
columns) 

i  (0) 

Initial  conditions 
on  the  state  vector 

8F1G.6 

(read  from  top  of 
column  to  bottom) 

i 

a 

! 

Table  II-l:  Definition  and  Format  of  Data  Cards  for  Program 
0PC0N1. 


START 


oil 


Figure  II-2:  Flow  Chart  Of  Program  OPCON1 
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al 


FIGURE  II-2:  (continued) 


1* 


a2 


_ 3[ _ . 

PRINT: 

P,iJ),and(Q+P) 


SM(I, J,k)-  Z(M-I, J ,k) 


W 

PRINT: 

INVERTED  (Q+P) 


a3 


Figure  II-2;  (continued) 


15 


a3 


STOP 


Figure  II-2:  (concluded) 
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START 


Figure 


II-3:  Flow  Chart  For  Subroutine  PROD (A, B,C,L>M, N) 


START 


Figure  11*4*  Flow  Chart  For  Subroutine  SUM(A>B,CsN) 


START 


Figure  II-5  :  Flow  Chart  For  Subroutine  TRANSQ(A, B,N) 
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START 


Figure  II-6I  Flow  Chart  For  Subroutine  TRANCOL(A, B,N) 


START 


Figure  II-7:  Flow  Chart  For  Subroutine  Money (MON, E,Q,N) 
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EXIT 


Figure  IX-8:  Flow  Chart  For  Subroutine  PPSI (B ,C ,D»Ej F >N) 


2  1 


EXIT 


Figure  II-9  :  Flow  Chart  For  Subroutine  PP(AjB>C,N) 


EXIT 


Figure  11*10:  Flow  Chart  For  Subroutine  FFtA^B^CjDjN) 


CHAPTER  III 


PROGRAM  OP CON 2 

Program  OPCON2  was  written  as  the  final  phase  for  a  versatile  simu¬ 
lator  to  determine  a  constant  gain  matrix  for  optimization  of  a  control 
system  according  to  the  following  cost  function: 

J (N)=minimum  ?  (x7(k)Qx (k)+ru^ (k-1 ) )  (45) 

k=l 

Again  u(k)  is  a  scalar  control  and  the  Q  matrix  is  determined  in 
the  same  manner  as  outlined  in  the  previous  chapter.  The  variable,  r, 
is  a  weighting  cons  tan t , (Lagrange  multiplier)  determined  by  user  needs 
and  acts  as  a  quadratic  integral  constraint  on  the  control  effort;  such 
as  a  limit  on  energy  necessary  to  provide  the  desired  control.  By  let* 
ting  r  equal  zero,  we  have  the  cost  function  of  0PC0N1  and  identical 
results  are  obtained  from  the  two  programs.  With  the  use  of  0PC0N1 ,  no 
restraints  are  put  on  the  control  required  by  the  system  which  may  re¬ 
sult  in  the  necessity  of  an  unlimited  supply  of  fuel  to  obtain  the  de¬ 
sired  optimal  results. 

Again  let  us  consider  the  system: 

x  =F x  +  Du  (1T) 

whose  discrete  solution  is: 

,x(k+l)-  <j> x.  ( k )  + Au  ( k )  (4 r ) 

where  u(k)  and  A_  are  defined  in  the  same  manner  as  in  the  previous  chapter. 

By  starting  with  the  last  stage  first,  equation  45  becomes: 

J (N)=minimum  |  xT (N)Qx(N)+ru^ (N-l)j  +J (N-l)  (46) 

or  by  substitution  of  equation  4'; 

J(N)=minimum[($x(N-l)+Au(N-l))TQ(0x(N-l)+Au(N-l))+ru2(N-l)]  +J(N-1) 

1  J  (47) 


2  3 


We  now  minimize  J(N)  with  respect  to  u(N-l)  noting  that  J(N-l)  is  not  a 
function  of  u(N-l). 


3  J(N) 

3u(N-l) 


0  =  (N- 1  )+A^QAu  (N- 1 )  +ru  (N-l ) 


(48) 


or. 


u(N-l)=  -  ATQ» 

ATQA+r 


x(N-l) 


(49) 


As  in  Chapter  II,  the  denominator  of  equation  49  will  always  be  a  scalar 
number;  therefore,  no  matrix  inversion  is  required. 

By  proceeding  back  one  more  stage,  the  control  policy  u(N-2)  may  be 
optimized.  Equation  45  becomes: 

J (N)=minimum  (N) Qx (N)+jcT  (N-l)Qx(N-l)+ru^  (N-  1)  iru^  (N-2)^J+J(N-2)  (50) 

The  following  symbology  will  be  used  for  clarity  in  this  minimization: 

A  4  4>x(N-2)+Au(N-2)  (51) 

and, 

A TQ<f) 


B  & 

T 

A  QA+r 

By  direct  substitution  using  4f,  51,  and  52,  equation  50  becomes: 

J  {  N)  =minimum  {  ( <j>A+ABA  )  TQ  ( (j>A+ABA)+ATQA+r  (BA) T  (BA )  +r  u  2  (N-  2  )  j  +J  (N-  2  ) 
By  rearranging  terras, 

J(N)=ininimutn|ATQA+((,<t>+AB)A)TQ(('<))+AB)A)+r(BA)2+ru2(N-2)j  +J(N-2) 

We  now  let, 

T  _  _ 

al  "  B  "  - 


(52) 


AtQA_+i 


(53) 


(54) 


(55) 


and 


Equation  54  now  becomes: 


^1=  <H‘Aa| 


(56) 


J  (N)=minimum^ATQA+AT^Qi{)iA+r(a^A)2+ru2  (N-2)j  +J(N-2)  (57) 


We  now  minimize  J(N)  with  respect  to  u(N-2) , 


9  J  (N)  t  t  T  T 

=0=  A  QA+A  ^Q^jA+r  a1Aa1A+ru(N-2) 


Solving  for  u(N-2)  by  substitution  of  51  into  equation  58, 

,T  T> 


u  (N-2)  =  -  /ATQ+AT^Qi})i+rATa1ai  | 

[  ATQA+AT^T1Q^1  A+r  ATa  l  A+r] 


x(N-2) 


We  now  let. 


P1=<f’iQ<l'i+Q+rala1 


Equation  59  becomes: 


u (N-2)  =  -  ATpi<l> 


ATpiA+r 


x (N-2) 


As  we  did  in  the  previous  iteration,  we  let, 


a  l  =  -  % 


A^p  A+r 
~  1“ 


and , 


(58) 


(59) 


(60) 


(61) 


(62) 


*2-  <HAa2  (63) 

By  proceeding  back  one  more  stage,  a  recursive  relationship  may  be  defined. 
We  shall  define  the  following  quantity  for  brevity  in  this  development: 


G  £  (j>x(N-3)+Au  (N-3) 


(64) 


The  cost  function,  45,  may  be  written  as, 

J (N) =minimum (N)Qx(N)+x^ (N-l) Qx(N-l)+x  (N-2)Qx(N-2) 

+ru2(N-l)+ru2(N-2)+ru2(N-3)J+J(N-3)  (65) 

By  substitution  of  4',  55,  56,  62,  63,  and  64,  equation  65  becomes: 
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J  (N) “minimum  j  (<^2 ^ iG) TQ  (^2*  1G ) +r  ( ^ 2^ i  £G)T  <^2alTG) 

+  (^2^)^('('2^)+r Ca2  G)^(a2  ^) 

+GTQG+ru2(N-3)j  +J(N-3)  (66) 

Minimizing  the  cost  over  N  stages  with  respect  to  u(N-3), 

3  J(N)  „  .T.T.T  .  .  ,T , T  .  ,T 
3u(N-3)  =  0  =  -  QG 

+rAT|^2aia^2^+r^ra2a2^+rU^~^^  (67) 

Substituting  64  into  equation  67  and  solving  for  u(N-3),  we  find  the 
optimum  control  policy: 


u(N-3)  = 

{  AT^lQ^1^2+— T^2^^2+— TQ+r— T*2ala1^2+r-Ta2a2  )  ^ 

^AT^2^^i^2~+~T^2^2“"*"~T^~+r“T^2aiai^2~+r~Ta2a2~+r]  ^ ^ 


Combining  terms, 

u(N-3)=  -  { AT^2(^TQ^l+Q+raiaI)|/'2+— TQ+rATa2aI  1  ^ 

(^lQ^l+Q+ral^I)  li^A+A*’ QA+rAT a2a2A+r] 


x(N-3) 


(69) 


From  equation  60,  equation  69  becomes: 

new-'*'!-  -  f  Ar^2Pl’J'2+ATQ+rATa„a„T  ? 

u(M-3)-  -  t-  2  w  ~  y  ~  ?  2  J -  x(N-3)  (70) 

{  AT^2pi  ^2A+ATQA+rATa2  a2T-+r } 

We  now  letj 

P2=^2PL'|/2+Q+ra2a2T  C71> 

and  equation  70  becomes: 

u(N-3)  =  -  r~Tp2^  x(N-3)  ‘  (72) 

A  P24+r 

Let  j 
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(73) 


T=  -  AV* 


ATP  A+r 


and , 


+4Aa3 


(74) 


As  in  Chapter  II,  a  recursive  relationship  may  now  be  shown  without 
proceeding  back  another  stage  and  may  be  shown  to  be  true  by  an  induc¬ 
tive  proof*  Let 


and , 


where , 


aT  =  -  -  k-l$ 

k  ATpk-l-+r  '  aT  -  0 


t  =  (|)+&a  ,  <|>  =  0 

k  “  k  J  o 


P,  =  ^Pi  +Q+ra  a  T,  P  =Q 
k  k  k-1  k  k  k  0 


(75) 


(76) 


(77) 


and , 


u(N-k)"  a  x(N-k) 
k“ 


(78) 


from, 

2c(kfl)=  $x.(k)+Au(k)  (79) 

By  letting  r  equal  zero  and  solving  the  same  illustrative  example  as 
in  the  previous  chapter,  the  same  results  are  obtained  as  illustrated 
below* 


and , 


and , 


a  T =  JLAS*.  r_2.0  -2.0 

1  T 

&_  QA+r  L 


iJj  =  0+Aa, 

1  “  1 


0  0 

-2.0  -1.0 


Pj=  ^Qi|/1+0+ra1a1T=Q= 


1.0  0 

0  0 


(80) 


(81) 


(82) 


2  7 


(83) 


and  therefore, 

u(k)=a^T>c(k)=  -  2x^  (k)-2x2  (k) 

Equation  83  indicates  the  same  results  as  in  Chapter  II  and  perscribes 
the  same  phase  plane  (Figure  II-l). 

Through  the  use  of  the  recursion  formulas  75,  76,  77,  78,  and  79, 
PROGRAM  0PC0N2  was  written.  Figures  III-l  through  III-5  and  Figures  II-3 
through  II-6  represent  abbreviated  flow  charts  of  the  program  and  the  sub¬ 
routines  used  in  the  program. 
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START 


al 


Figure  III-l :  Flow  Chart  for  PROGRAM  OPCON2. 
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al 


a2 


Figure  III-l:  (continued) 
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Figure  III-l:  (continued) 
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*  f 


FIND: 

u(l,l) 

z=  u(l,l) 

L  =  1 

] 

i 

FIND: 

DOL 

TCOST  =  TCOST+DOL 


V 

a4 


Figure  III-l:  (continued) 
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a4 


_ 1 

l _ 

PRINT: 

L,Z,X(1),  and  TCOST 

1 

_ 

KN=  1 

Y(J, 1)=  X(J,1) 


'  j 

N  ' 

< 

1 

>/ 

t 

X=<f>*Y 

■ 

DELTA (J , 1 )=Z*DEL ( J , 1 ) 
X(J,1)=  X(J,1)+DELTA(J,1) 
AT(1,J)-ALFAT(I,J) 


J  :N 


c*5 


Figure  III-l: 


3  3 


(continued) 


22 

A 


a6 


A 

22 


Figure  III-l:  (continued) 
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STOP 


Figure  III-l:  (continued) 


3  5 


START 


EXIT 


Figure  III-2:  Flow  Chart  for  SUBROUTINE  A TRAN (AT , P , PHI , DEL ,  R ,  N) 


START 


Figure  III-3:  Flow  Chart  for  SUBROUTINE  PPSI(PSI,PHI,DEL,AT,N) 


START 

Y 


FIND: 

YT 


A  =  YT*Q 
AA  =  A*Y 

D0L=  AA(1 

,1)+R*Z*Z 

A 


EXIT 


Figure  III-4:  Flow  Chart  for  SUBROUTINE  COST(DOL,Y,Q,R,Z,N) 
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START 


o?  0 

Figure  III-5:  Flow  Chart  for  SUBROUTINE  PP(P,PSI,P1,Q,AT,R,N) 
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>/ 


AE=  AC+Q 
P=  AE+AD 


A 

EXIT 


Figure  III-5:  (concluded) 
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Figure  III-l  shows  that  the  inputs  to  the  program  must  be  read  from 
data  cards.  Table  III-l  shows  the  order  in  which  these  cards  are  read, 
the  format  required  on  each  card,  and  the  definitions  of  the  variables 
concerned. 


Variable 

Definition 

Format  of 

Data  Card 

M 

Number  of  stages 

Plus  one 

19 

N 

Order  of  system 

19 

IT 

Frequency  of  the 
desired  printing 
of  output 

19 

r 

Weighting  constant 

of  cost  function 

1 

F10.6 

Column  vector  of 
system  impulse 
responses  due  to 
control  inputs 

8F10.6 

(read  from  top  to 
bottom  of  column) 

♦ 

State  transition 
matrix 

nF10. 6 

(read  in  a  row  at 
a  time) 

Q 

Weighting  matrix  of 
the  cost  function 

nF10 .6 

(read  in  a  row  at 
a  time) 

x(0) 

, _ 

Initial  conditions 
on  the  state  vector 

8F10 , 6 

(read  from  top  to 
bottom  of  column) 

i  i 

Table  III-l:  Definitions  and  Format  of  Data  Cards  for 

PROGRAM  OPCON2 . 


t 
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The  output  of  PROGRAM  0PC0N2  is  the  value  of  each  state  variable, 
the  control  signal,  and  the  total  cost  of  all  previous  stages  includ¬ 
ing  the  present  stage.  The  amount  of  printed  output  may  be  varied  by 
the  use  of  the  IT  variable's  input  data  card.  If  all  output  is  to  be 
printed,  then  IT  must  be  set  to  one  by  the  data  card.  If  the  user 
wishes  to  have  1,000  stages  and  has  a  requirement  of  print  out  of 
only  a  sampling  of  these  stages,  IT  should  be  set  accordingly.  By  let¬ 
ting  IT  equal  ten,  only  every  tenth  stage  will  be  printed  out.  Another 

T 

valuable  output  is  the  value  of  a  which  approaches  the  constant  gain 
matrix  necessary  for  optimal  control  of  a  continuous  system.  This  is 
achieved  by  using  a  small  sample  period  and  many  iterations  thereby 
approaching  the  simulation  of  a  continuous  system-  0FC0N2  is  limited 
for  this  use  because  the  computer  memory  size  limits  the  program  to  ap¬ 
proximately  2,000  stages*  If  the  user  has  no  need  for  the  outputs  of 
the  values  of  the  control  signals,  state  variables  and  total  costs, 

then  the  program  is  unlimited  in  its  number  of  stages  because  the  value 
T 

of  a  may  be  printed  out  as  often  as  desired  without  storing  it  in  mem- 
K 

ory.  This  alteration  has  been  made  in  PROGRAM  GPC0N3,  and  a  listing 
of  this  program  may  be  found  in  the  appendix. 

The  versatility  of  this  program  may  be  further  shown  in  using  it 
to  simulate  an  optimum  control  in  accordance  with  the  cost  function, 


J  (N)=  minimum  |^XT  (N)X(N)j 


(84) 


which  we  shall  call  the  final  value  criterion.  Using  the  same  method 
of  derivation  as  before,  we  find  the  following  recursion  formulas: 
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ip  *<H-AaL  ip  —  0 

K.  IS.  O 


(85) 


and , 


K 


T 

*  PM* 

— T^k-  1— 


T 

a 

o 


=0 


where , 


P  =  a1  p  it  ,P  =  I 
K-l  \-l  K-2  K-l  ° 


and. 


u(N-K)=  a  a(N-K) 
K 


(86) 


(87) 


(88) 


To  use  0PC0N2  for  this  cost  function,  we  need  only  to  set  P  equal 
to  the  identity  matrix,  Q  to  the  null  matrix,  and  r  equal  to  zero. 

As  an  example  of  this  cost  function,  we  shall  again  use  the  same 
illustrative  system.  We  find  that: 


and , 


and , 


and , 


T 


.4 


-.4 


pi  =  +r 


.8 

,4 


.4 


.2 


T 

a 

2 


ATP1<j> 

ATP1fi 


(89) 


(90) 


(91) 


(92) 


therefore , 


u(0)  =  -  X1(0)  -  1.5X2(0) 


(93) 
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Given  initial  conditions  of  X^(0)  =  1.0  and  X2(0)  =  0.0,  equation 
93  results  in, 

u(0)  =  -1.0  (94) 


therefore , 

Xx(l)  =  X1(0)+0,5u(0)  =  0.5  (95) 

and, 

X2(l)  =  X2(0)+u(0)  =  -  1.0  (96) 

From  equation  89, 

u(l)  -  -.4X1(1)  -  1.2X2(1)  =  1.0  (97) 

and , 

* 

Xx(2)  =  X1(l)+X2(l)+.5u(l)  =  0.0  (98) 

and, 

X2(2)  =  X2(l)+u(l)  =  0.0  (99) 


Figure  III-6  is  the  phase  plane  resulting  from  the  calculated 

J 

controls  and  shows  that  the  final  states  at  N  equal  two  are  definitely 
at  a  minimum  and  the  cost  function  equals  zero. 
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With  the  use  of  0PC0K2  and  GPC0N3,  many  control  systems  may  be 
simulated.  Using  these  two  programs,  a  system  may  be  designed  for  a 
number  of  optimization  criteria  while  greatly  reducing  design  time. 
In  the  next  chapter,  a  typical  system  will  be  considered  assuming 
an  optimization  criterion  and  the  resulting  design  will  be  obtained. 
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CHAPTER  IV 


EXAMPLE  PROBLEM 

To  show  the  use  of  the  programs  that  were  written  and  developed 
in  the  preceding  chapters,  the  following  example  problem  was  chosen. 
The  transfer  function  is  typical  of  a  servo  motor,  the  roll  attitude 
control  of  a  missile,  and  many  others.  Figure  IV-1  shows  the  open 
loop  plant  to  be  controlled. 


Figure  IV-1:  Block  Diagram  of  the  Open  Loop  Control  System. 
The  following  state  equation  defines  the  system: 


0  1 

0 

= 

X  + 

£  ~L 

l 

j 

(100) 


With  the  use  of  PROGRAM  PHIDEL  and  a  sample  period  of  one  second, 
and  A  are: 


1.00 

.631 

.361 

x(k+l)  = 

_x(k)  + 

it-0 

.  368_ 

.632 

“00 


(101) 


To  show  the  flexibility  of  PROGRAM  0PC0N2,  four  different  cost 
functions  will  be  considered  and  compared.  The  first  cost  function 
will  be  the  terminal  value  function  for  time  optimal  control.  (Case  I) 


T 

J(N)  =  minimum  xA(N)x(N) 


(8V) 


This  will  then  be  compared  to  the  terminal  value  control  with 
an  energy  control  shown  by  equation  (102).  (Case  II) 
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(102) 


J(N)  =  minimum  x^(N)x(N)  +  ^~  u^(K-l) 

k=l 

The  next  case  of  interest  minimizes  the  sum  of  the  squares  of  the 
states  with  an  energy  limitation.  This  is  represented  by  equation  (45'). 
(Case  III) 

k=N 

J  (N)  =  minimum  acT  ( K)  Qst  (K)  H-trii  ^  (K— 1 )  (45') 

k=l 

where , 

Q  =  I  and  r  -  1.0  (103) 

The  final  cost  function  will  also  be  that  of  minimizing  the  states 

with  an  energy  limitation;  however,  by  use  of  PROGRAM  0PC0N3,  a  steady 
state  feedback  matrix  was  found.  These  values  were  then  used  in 
PROGRAM  0PC0N2  to  obtain  a  transient  solution.  (Case  IV)  The  elements 
of  this  control  are: 

aT(l,l)  =  -0.627797  (104) 

and , 

aT(l , 2)  =  -0.615037  (105) 

With  the  use  of  PROGRAM  GPCGN2*  a  transient  solution  was  obtained 
in  each  of  the  four  cases*  Case  II  was  computed  for  a  twenty  second 
trajectory  in  the  phase  plane s  Figure  IV-3*  Cases  III  and  IV  were  each 
run  for  ten  second  trajectories  while  Case  I  was  optimized  in  two  samples* 
Figure  IV-2  is  a  plot  of  control  versus  time  for  each  of  the  cases* 

It  is  of  interest  to  note  that  both  cases  III  and  IV  utilize  the  same 

controls  during  the  first  seventy  percent  of  their  trajectories;  only 

when  the  control  is  very  small  in  the  last  couple  of  sampling  points 
does  the  control  required  differ*  This  thereby  leads  to  identical  trajec¬ 
tories  in  the  phase  plane  indicating  that  constant  gain  amplifiers  in  the 
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feedback  loop  are  just  as  effective  as  variable  gain  amplifiers  for  com¬ 
pensation  of  this  plant.  Case  II  requires  nearly  a  constant  control  ef¬ 
fort  after  the  initial  control  jump  has  been  made.  The  required  control 
for  the  twenty  second  trajector  is  approximately  equal  to  one  half  of 
the  necessary  control  for  a  ten  second  trajectory. 

Figure  IV- 3  is  a  phase  plane  plot  of  the  four  cases.  One  thing  to 
note  here  is  the  trajectory  of  case  II.  As  the  total  number  of  sample 
periods  is  increased*  the  necessary  control  is  decreased  as  indicated 
above.  This  results  in  the  trajectory  approaching  the  axis  as  the 
number  of  sample  points  increases  indefinitely.  This  means  that  the 
cost  function  used  is  optimized  when  the  final  value  of  the  states  is 
zero  and  the  control  is  nearly  zero.  The  point  in  the  phase  plane  would 
drift  in  requiring  only  a  very  slight  amount  of  control  to  stop  it  at  its 
steady  state  value. 
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Figure  IV-2:  Plot  of  Control  Versus  Time  for  the  Four  Cases. 
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Figure  IV-3:  Phase  Plane  Plot  of  the  Four  Cases. 
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APPENDIX  I 


PROG 7  A*';  PH  I  DEL 


PROGRAM  PH  I  DEL 

QD I  MENS  I  ON  F ( 1 2  » 1 2 ) *  PH  I ( 1 2  * 12  >  ,  TERM (12*12) * WORM (12.12)* 
1  DEL (12) »DELK( 12,12) » TELiM  ( 12  *  12  )  ,DELP(12»12) ,D(12) 


C 

r 

r 


c 

0 


THIS  PROGRAM  CALCULATES  THE  STATE  TRANSITION  MATRIX*  PHI-  AND  1  HE 
EXTERNAL  FORCING  FUNCTIONS  IMPULSE  RESPONSE  MATRIX,  DEL*  FROM  THE 
“'•TRICE'S*  F  AND  D,  THIS  IS  ACCOMPLISHED  ?Y  USING  THE  TAYLOR  S'  RIDS 
EXPANSION  FOR  THE  EXPONENTIAL  FUNCTION.  THE  CALCULATION  CDR  THE 
DEL  f.TRIX  IS  GOOD  FOP  ONLY  "ON  TIME  VARYING  DYNAMICS  RECAUSE  THE 
INTEGRATION  IS  ACCOMPLISHED  WITHIN  THE  PRQORA"  FOR  THIS  CASE-  ONLY. 

f 

PHI  =EXP*«  (  F*T  )  -  1  -M  F*T  )  +  (  (F*T)**2)/I  1  •  "  ■  2  •  P  )  4 - +  (  IF*T)*«.N)/(  1.0*.  .  *N  ) 

DEL= INTEGRAL  FROM  ZERO  TO  DT  OF  PHI { DT- T AU } *D 2 S T A U 
THE  REQUIRED  INPUT  DATA  CARDS  ARE, 


(A)  N  (SIZE  OF  THE  INPUT  MATRIX.  FORMAT  IE) 

(B)  DT  (SAMPLING  PERIOD,  FORMAT  1  FI  5. 12) 

(C)  TEST  (TEST  FOR  VAX  I  HUM  ALLOWABLE  ML'MjjEP  I M  A  TERM  TO  TRUNCATE 

THE  SERIES,  'rOR“M  1F1S.12) 

ID)  F  (DISTRIBUTION  MATRIX,  FORMAT  MUST  -c  INSERTED*  STATEMENT 
ML1 1 '  R  -  rj  3 1  ,  *■' r  1 n  •  6  1 

(E)  D  (COLUMN  VECTOR  RELATING  THE  CONTROL  VARIABLES  TO  THE  5YSTE" 
DYNAMICS,  ^0R*'4T  S“ln.5) 


OTHER  FORMAT  CARDS 


WHICH  MUST  3E 


INSERTED  ARE  MRS  34  AND  35  WHICH 
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PROGRAM  PHI  DEL  (CONTINUED) 
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APPENDIX  I  i—2  PP  OG  ?  £  r  P u I D  E  L 


Ln 

U> 


DC  50  10=1} 
DO  5A  I  0=1} 
50  JM=1 , 


r\ 


4  1 


DELHI IR , 

IC) 

- 

”  Lj‘ 

EL1-'! 

!  S, 

IC  )- 

tel: 

M  ( 

I  R  } 

JN 

WORM ( I 

O 

IC  ) 

=TERK( 

I  R. 

JN  )  "•  F  [  J 

■i . 

IC) 

.H.  r 

DO  41 

IP 

=  1 } 

M 

DC  41 

!  C 

=  1  » 

\j 

TERM ( I R » 

IC  ) 

-  j 

-  P  M  ( 

IP} 

IC) 

TELM ! I 

R  } 

I  C  ) 

EL’!  ( 

I  R } 

IC  ) 

DELP ( I 

P  } 

IC  ! 

=  r'£LP;' 

IP} 

I  C  )+TEL 

M  ( 

i  R } 

I  C 

DFLMI I 

R  > 

IC  ) 

=  0 

.0 

WORK ( I 

R  t 

IC  ) 

-o 

.  3 

M  =  TM 

PRINT 

34 

)  M  • 

f  ( 

7  e  f?  y 

(  IR 

;  IC) 

}  I  c 

"1 

}  N  ) 

:  I 

( JN»IC)*DT/(  T  ’!■ 
<-  i;;0  RM  {  I  R  *  T  C  ) 


1.0) 


SUKATI ON  OF  T  ME  7ER*'S  IN  THE  SERIES 


DO  51  I R= 1 } u 
DO  51  IC=1,N 

51  PHI ( IR,  I  C  )  =  PH  I  (  IR>  1  C  )  +  T_  RK  l  I  R  •  I  C  ) 


C  TEST  FOR  THE  TERMINAL  T  rR“  IN  THE  SERIES  A 'O’  PRINT  T  HE  PHI  MATRIX 


ABC=0. 0 
DC  42  IR~1,N 
DO  42  IC=1}N 
AA=T£RK ( ! R }  I  0 ) 
45=425-  (  V\) 

IF  (  ABC-Ab)  43.*42  *42 
43  ASC=A3 
42  CONTINUE 


r 

4 


r i 

«/ 


APPENDIX  I 


dPOSp:>-'-'  OH  I  DEL  (  CONT  I  N'JE:~> ) 


IF  (  TEST- ADC  )45  ,45 ,46 

45  GO  TO  44  ' 

46  PRINT  35  ( ( PHI ( I R * IC ) » I C=1 *N ) » I R=1 . N > 


CALCULATE  AND  PRINT  THE  DrfL  MATRIX 


DO  5  2  1=1  ,  N 
DO  5  2  K- 1 » N 
DO  52  J  =  1 ,  N 

5  2  DEL!  I  )  =DEL  (  I  )  +PHI  (  1  »  J )  *DELP  (  J  ,  Kl  *  D  <  K ) 
PR  I  NT  5  3  ! DEL ( I ) » I =1 »N ) 

31  FORMAT  (8F1Q.S) 

32  FORMAT  (15) 

33  FORMAT  ( 1F15.12 ) 

34  FORMAT  ( / // 7HT T ( I , J ) , I  2/ ( 3F1 2 . S )  ) 

3  5  FORMAT  ( / / / /OX , 8HRH I  (  I* J  )  / / / ( 3  FI  2 . 8  )  ) 
5  3  FORMAT  ( /  /  / /9X , SHDSL ( I) / / ( 8F 1 5 . ? ) / / > 

5  4  FORMAT  ( / / , 7HF (  I  ,  J  !  = / /  .  !  3 F 1 0 . 4 )  J 
5  5  FORMAT  (  /  /  »  5HD( I ) =// * ( 8F 10.4 ) ) 

56  FORMAT  (// . 3HDT=. t 1F13.3 ) ) 

END 

END 


APPENDIX  I 


I  -4 


PROGRAM  PHIDEL  (CONCLUDED) 


u  u  O  U  W  u  u  u  u  u 


APPEND I  X  II 


PROGRA'-1  0PC0N1 


PPOG^AV'  OPCOM1 

OD I  MENS I  ON  A ( 8*8 } ,3 ( 3 , 8 ) »Q( 3*8 1 ,D(8 ,8 ) >X ( 8 ,40 )  ,G(40) ,H<  8  »8l , 
1P<40» 8*8 ) ,R(40*3,S) »Z (40,8,3) »  COST  t  40 ) ,Y(8,8)  *ZZ(S,3)  * 

2RR f  8  >8 )  ,  P 1  (  6  ,  S  ) »  S  M  (  4  0  •  8  »  8  ) *  F  (  4  0  •  3  *  8  ) »  A  A  (  8  *  8  )  » XX (8*8)  »E ( 8  ,  8 )  * 

3 MON ( 8  » 8 ) » SMM t  8  *  8 ) , U ( 40 ) , UU I  8  » 3 ) 


THIS  PROGRAM  IS  DERIVED  USING  THE  COST  FUNCTION,  J ( N ) =M I N I  MUM ( SUM 
X (N >T*Q#X ( M) ) .  THE  PROGRAM  WILL  YIELD  OPTIMUM  CONTROL  FOR  UP  TO  AN 
EIGHTH  ORDER  SYSTEM  OVER  40  STAGES.  THE  INPUT  DATA  CARDS  ARE, 

(A)  M  ( NUMBER  OF  STAGES  PLUS  ON-,  FOR'' AT  15} 

(SI  N  (ORDER  Oc  SYSTEM,  FORMAT  15) 

(C)  DEL  (COLUMN  MATRIX  READ  IN  ORDER  BY  ROW,  FORMAT  8 FIG. 6) 

(D)  PHI  (TRANSITION  MATRIX  READ  IN  A  COLUMN  AT  A  TIME,  FORMAT  8F10.6) 

(E)  Q  (NXN  WEIGHTING  MATRIX  READ  IN  SAME  AS  PHI,  FORMAT  8F10.6) 

IF)  X  (COLUMN  VECTOR  OF  INITIAL  STATE  CONDITIONS,  FORMAT  BF10.6) 

THE  FOLLOWING  RECURSIVE  EQUATIONS  WERE  DERIVED  AND  ARE  USED  AS  THE 


BASIS  FOR  THE  OPTIMIZATION  OF  CONTROL, 

C  PSI  IK+1 )  =  ( I - { DEL * DELT  * { 0+  P ( K )  ))  /  ( DELT  * <  Q  +  P  f K )  ) *  DEL ) )*PHI  I  1 ) 

C  P (<-l ) =PSIT( K-l )*P( K-2 H PSI ( <-l  )  (2) 

C  U  (N-K  >  =-(  (  DELT*(G+Pt,K-l  )  )  *PHI  )  /  1  DELT*  I  Q+P  I  K~1  )  )  *DEL  )  )#X(N-K  )  (  3  } 

C  X { K+ 1 J  =  PH I *X ( <  )+DEL*U ( K )  (4) 


APPENDIX  II 


I  I-I 


PROGRAM  CP CONI  (CONTINUED) 


C  E  O'.iAT  IO^S  1  A D  2  CONSTITUTE  SUBROUTINES  IN  THIS  P^OGRA'-’,  ANOTHER 
C  SUBROUTINE  GENERATES  THE  RESULTING  COEFICIENT  MATRIX  FOR  X ( N - < )  IN 
C  EQUATION  3.  ( SUBROUT I NF  Fpj 


ui 

O' 


71 


72 


READ  202 
READ  202 
READ  201 
READ  201 
READ  2 C 1 
READ  201 


P’  ) 

( N  ) 

(A(Ij])sI=l;N) 

{  I  B(  I  *  J)  )  I  =  1  -  N  I  •  J-  1  • ' '  1* 
t  t  Q  {  I  5  J  }  ■  I  =  1  *  N  )  *  J  =  1  j  '•: ) 

(  X  [  I  .  1  )  ,  I  =  1  ,  N  ) 


PRINT  2  02  (  ) 

PR  1 NT  2  02  (  N  ) 


PRINT  201  (  M  I ; 1 )  *  I  =  1 j  N ) 

PRINT  2C1  (  (3(1 •  J  ) *  I  =  1  •  N  j  j  J  -  1 .  N  I 

PRINT  201  (  (Q(  I  » J)  j  !=1  *.M)  f  J=liN) 
PR  I  MT  2^1  (  X  (  I  .  1  )  *  1  - 1  -  N  I 
CALL  DELTA  {  D  .  A  >  N ) 

DQ  71  1=1. N 

DC  71  J=1 »v 
2(1)1. J)-Q(I;J) 

ZZ(  I  * J 1=2(1 , J  . J) 

DO  ^  2  I  =  1  P'  - 1 

CALL  P.PS  I  (  D  )  ZZ  :  A  •  7  •  R  p  *  N  ) 

CALL  dd  (  ° D  *  Q ■■  c'  1  .  N  ) 

CALL  SU"( 7 2 ■ D - PI ■ N } 

DD  72  J=1,M 
DO  72  K- 1 ) N 
P(  I  .  J.<)=P1  (  J.O 


R ( I  * J.K ) -RR ( Jj  K ) 

Z( I X1 » J)<)=ZZ  t  J)X ) 

PRINT  10Q4 

PRINT  K  2  (  (  <  I  .  J »  <  ?  R  (  I  *  J  )  K  }  j  .<=  1  *  N  )  *  J  =  1  ■  N  )  >  I  =  1 » M  ) 
PRINT  lrn5 

P?  I  NT  1 "  2  (  (  (  I  »  J  ,  '<  -  P  (  I  .  J  .  <)  •  :<=  1  ,  N  )  ,  J=  1  ,  N  )  .  I  -  1  .  M  ) 


APPENDIX  II 


I  1-2 


PROGR  A  v  0°  C  0  N 1  (CONTINUED) 


L_n 


PRINT  1006 

PRINT  102  ( t ( I > J.K»Z (I , J , K=I ■ M) , J=1 , N) , 1=1 .  Y ) 
DO  73  1=1 
DO  73  J  =  1  ,« 

DO  73  K  =  1 » N 

73  S  ‘  “  <  I  •  J  *  K  )  =  Z ( N-I , J; <) 

PRINT  1007 

PRINT  102  (  (  (  I  »  J  •  K  •  S1-'  (  I  f  J  *  <  )  »  <=  1  ■  N  }  ■  J=  1  5  N  )  t  I  =1  t  M  ) 
PRINT  1002 
DO  74 

DO  75  J  =  1 » N 
DO  7  5  K  =  1 >  N 
7  5  S M  (  J  » '<  )  =  Sv  {  I  »  J  *  K  ) 

CALL . Fr ( H  *  A  »  ? i SMM  ■  N ) 

DO  2  0  J  =  1 , N 
Rn  2 0  <  = 1  ,  N 
2  0  H(  J,<] 

no  75  J  -  1  ?  M 

DO'  75  -  1  , M 

7  6  Y  (  J  t  K )  =  0 .  P 
DO  79  J=1,N 
79  Y ! J  ,  1 1 =X ! J,  I  ) 

CALL  PROD  t  UU • H , Y  *  1 ! 1  - N ) 

U (I) =UU (1,1) 

C ^  L  L  D ROD ( *  A , 4 ! UU • N •  1 ,  1  ) 

CALL  PROD  (  XX  ■  3  *  Y  O)*  1  ,  1 

DO  77  J  =  1  ,  N 

7  7  XI  J,  1  +  1  )  =XX  (  J;  1  )  +  AM  J,  1  J 
DO  7  8  J  =  I » N 
F  (  I  ,  1  ,  J  )  =  H  I  1  ,  J  ) 

78  E(J»l)=X(JsI) 

CA  L  L  *'ONr  Y  ( ‘ ' Q N  ,  -  ,  0  ■  N  ) 

PR  I  NT  100  (I  ,M0N( 1,11) 

7  4  COST  (  I  )  =-MQN  (1,1) 

U(M)=0.C 
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CALL  PROD  (  a  A  ,  G  ,  C  >  N  >  N  ,  M ) 

CALL  PROP  I 9B,3, AA*N ,N,N) 

CALL  TRANCOL  t  D*PPi  N) 

CALL  PROD  (  Cf»'‘»r')N*  1  *N) 

CALL  PROP  { EE»0D*CC» 1 » 1 »N) 

DO  60  1=1  ,  M 
DO  60  J=1 ,M 
60  A  (  I  *  J )  =  C  .  0 
DO  6  1  I  =  1  ,  K 
<61  A  (  I  *  I  J  =  1  •  0 
DO  62  1=1 »N 
DO  62  J=1.N 

62  BB  t I  *  J  )  =  F  B  (  I • J1  /f E ( 1  - 1 > 

DO  63  1=1,0 
DO  53  J  =  1,M 

'63  BB(I,J)=MI,J) -QB  (  I  *  J  ) 

CALL  PROP  (  F,B3  ,E  •  0,  v,  0) 

END 

C  THE  FOLLOWING  SUBROUTINE  CALCULATES  P. 

SUBROUTINE  P  P  (  A  *  3  *  C  *  N  } 

D I  HEMS  ION  A <  S  » 8 ) » 9 ( 8 , 8 )  - C ( B ■ 3 )  -  0  {  8 • 3 ) »  A  A ( 8 ■ 8 ) 

CALL  T  RAN SO  t  A , A  A • M ) 

CALL  PRODfD-  B*  \.N  »  N*.N  } 

CALL  PROP! C* AA, D,N- N »N ) 

END 

C  THE  FOLLOWING  SUBROUTINE  CALCULATES  THE  COEF.  OF  THE  STATE  VARIABLES 
C  FOR  THE  CONTROL  POLICY. 

SUBROUTINE  FF ( A , 3 . C , D , N ) 

OD I  MENS  ION  Af8*S)»3(8»8)*Ct8*8)»D(8*B)rE(B*9)»S3(8»8)  *  G  (  8  *  8  )  » 

1H ( 8  ,8  ) 

CALL  PROD ( E  » D  >  C  »  M • N i  U) 


APPI 


DIX  I 


I  I  ~E 


PRC' 


OPCOMl  (CONTINUED) 


CALL  TRANCOL<3,93»SO 
CALL  PROD  (  A  i  02.  Z*  1  jNi  M ) 

CALL  PROD  (  GO- 9- N,  1  .  N) 

CALL  PROD ( H  »  B3 » S , 1 , 1 , N  ) 

DO  64  r  =1 »N 

64  A( 1 » I  ! =A( 1 , I ) /H ( 1 ,1  1 

C7M 

L  I  I  U 

C  THIS  SUBROUTINE  CALCULATES  THE  C05T  AT  EACH  STAGS. 

SU3R0UT  !  NS  ''0MEY  {  MQN  .  E  *0  !  M  ) 

DIMENSION  "OIK  S  ,5)0(30'!  01(3*  3  )  *  Q(  S»  S  )  »  HO  (8,3) 

CALL  TRANCOUEOTO!) 

CALL  PROD  CUD  » Q  »  E  >  N  j  1  i  N  ) 

S  CALL  .  PROD  C'ON » ET  »  "0«  1,100 

END 

C  THIS  SUBROUTINE  TRANSPOSES  ANY  CCLUME  MATRIX  (  1 ,  ‘J ) 

S'JB  ROUT  T  NE  T p  -  N C 0 L  (  "  ■  n  ,  N  ) 

D I  YENS !  ON  M3fS)tB(3*B) 

DO  12  I “ 1 s  N 
12  B ( 1  ,  I  )  -  M  I  ,  1 ' 

END 

C  THIS  SUBROUTINE  TRANSPOSES  A.  .NY  S  0  U  A  RE  “  A  T  R  I  X  !  ,N  X  ,N  ) 

SUBROUTINE  T R A N c 0 ( A , B , M) 

D I  HEMS  I  ON  A ( S , 3 ) »  3 ( 3  f  8  ) 

DO  .11  I  =1 ,  N 
DO  11  J  =  1  *  N 

11 

END 

C  THIS  SUBROUTINE  CALCULATES  THE  PRODUCT  OF  ANY  TWO  MATRICES  AS  LONG 
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PROGRAM  OP CONI  ( CONI ! NU ED ) 


AS  THE  PRODUCT  1$  ALGEBRAICALLY  PPr'DE 
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APPENDIX  II  I  1-7  PROGRAM  OPCOM1  (CONCLUDED). 


n  n  rv  n  n  n  n  n  n  n  n  n  n  n  n 


APPENDIX  I  I  I 


PROGRAM  OP COM 2 


PROGRAM  OP COM? 

ODI MENS  ION  AT ( 6,3 ) » ALFAT (200G»3)  » Q I  6 , 3 )  , P ( 3 » 3 )  ,  P 1 ( 3 , 3  )  , DE L ( 3  » 8  I  * 
IPS  I (a, 8) .PHI ( 6,8 ) ,X ( 3 ,3  >  , Y ( 8 ,3 )  .DELTA! 5 ,8) ,0 ( 3i8 ) 


THIS  PROGRAM  HAS 
x  {  n  )  t*c*X(  ksum 


3 i l i «  DERIVED 
R * 1 J  (  N-l  )  *2  ) 


USING  A  COST  FUNCTION,  J  (  I; )  =  M  I  N I  MUM  (  3U' 
THE  DIMENSION  STATEMENTS  LIMIT  THE 


>  E  R 


SYSTEM 


PROGRAM  TO  AN  EIGHTH  3RD 
FOLLOW  PNG  ITERATIVE  EQUATIONS  WERE  DEVELOPED 
ulLLMaNS  PA  INC  ±  PL ^.G  Or  DYNAMIC  Pi,  UGA  A  MM  I  hG  » 


/ 1 TH  2C00  ITERATIVE  STAGES.  THE 
A i .  DERIVED  USING  R  •  l  i 


A  I  K  }  T  =  ~(DELT*P(.<“1)  *  P  ii  I  )  /  (DELT^Pt  R- 1  )  *DEL+R)  <  i  ) 

PSI ( K) =PHI+DcL*A( < )T»  P  S I <  0 ) - C  (2i 


P  (  K  )  =  PS  I  t  X  )  T  P  (  K-  i  )  •"  r  S  i  (  <  )  ■*  0+  R  w  A  t  K )  *>  A  (  <)  T  »  P  t  0  )  -Q  t  3  ) 

X  (  <  )=PnI  ”  X  (X-l  J+DEL-U  (  <~l  )  (  A  ) 


U(N-<>=A(K)T*X(N-X) 


(5) 


EQUATIONS  i»  2»  AMD  j  CONST IIoTl  G  u  □  A  0  U  T I N l 3  IN  THIS  PROGRAM* 

THE  DATA  CARDS  REQUIRED  ARE  LISTED  IN  ORDER, 

(A)  M  (NUMBER  OF  STAGES,  FORMAT  15 } 

(d)  M  (ORDER  OF  SYSTEM,  FORMAT  13) 

(C)  IT  (MULTIPLE  OF  RESULTS  Or  STAGES  TO  GE  P.RIhTcD,  FORMAT  15) 


APPENDIX  III 


I  I  I  - 1 


PROGRAM  OP C Of42  (CONTINUED) 


c 

c 

(  D) 

R 

c 

c 

(  E) 

DEL 

c 

c 

(  F) 

PHI 

c 

c 

(  G) 

0 

r- 

V 

(  H) 

X 

c 

r 


CONTROL.  FORMAT  1F10.6) 

(THE  SYSTEM  IMPULSE  RESPONSE  MATRIX.  FORMAT  3F10.6,  READ 
IN  BY  COLUMN) 

(THE  S  Y  S  T  EH  T  RAMS  I T I  on  HA  T  RT  X  *  FORMAT  NF 1 Q . 6  »  READ  IN  by 
ROWS  ) 

(WEIGHTING  MATRIX  USED  IN  THE  COST  FUNCTION,  FORMAT  NF10.6, 
READ  IN  BY  ROWS) 

(INITIAL  CONDITION'S  ON  THE  COLUMN  VECTOR  Or  STATES,  FORMAT 
3F10.6,  READ  IN  BY  COLUMN) 

FORMAT  STATEMENT  NR  8  WILL  HAVE  TO  BE  ALTERED  DEPENDING  UPON  THE 
ORDER  OF  THE  SYSTEM.  ITS  FORMAT  SHOULD  BE  NF10.6. 


C  READ  DATA  AND  INITIALIZE  VARIABLES* 

2  READ  2<M> 

READ  2  IN) 

READ  2  (IT) 

READ  3  ( R ) 

READ  1  (  DEL  l  I » 1  )  » I =1*N  ) 

READ  8  ( ( PHI ( I , J )  , J=  1  ,N ) » 1  =  1  ,N ) 

R  r  AD  8  (  ( 0 (  I  »j)  »  J  =  1  ,  N  )  ,  I  =  1  *  N  ) 

READ  i  (X(I»1>»I=1,N) 

PRINT  2  (M) 

PRINT  2  (N) 

PRINT  2  (IT) 

PRINT  3  ( R ) 

PR  I  NT  I  < DEL (  1,1)  » I  =  1  .  N ) 

PRINT  6  {  (Phi (  I ,JJ ,J  =  1,N)  ,I  =  i,N) 
PR  I  NT  8  (  (  0  (  I  .  J  )  .  J=  1  ,  N  )  »  I  -  1 .  N  } 
PRINT  1  ( a (  I  ,1)  ,  I  =  1 , N ) 

TCDST=0.0 
KN  =  0 


APPENDIX  III 


I  I  1-2 


PROGRAM  OP  COM2  ( CONT INULD) 


DO  2  0  I  =1  »N 
DO  20  J=1,N 
20  PI  {  I  . J ) =QtI » J ) 


C  CALCULATE  AND  PRINT  A(A)T,  PSI(n), 


DO  23  1=1) M— ] 

CALL  A  T  RAN ( AT ,  P 1 ,  P H I  ,DEL»R»N) 

DO  2  2  J  =  1 ) N 

2  2  ALFAtHM-I  »J  »=AT  (1  *  J) 

CALL  PPG  I ! P  S I ) P  H I »  u  E  L  >  A  T  » N  ) 

CALL  PP ( P .PS I .  P 1 > 0 , AT , R, N ) 

KN=  KN+ 1 

IF  (RN-IT) 17,18,1b 
18  PRINT  4»I  » ( AT ( 1 , J I  »J  =  l»i\i) 

PRINT  5  )  I  »  (  (  PS  I  <  J  » ,<  >  ,  X=  1 » N  )  ,  J  =  i  , 
P  R I N T  6  » I ) (  (  P  (  J  » <  ) ,  A  =  1 » D )  ) J  —  1 )  j  s  ) 

KN  =  0 

17  CONTI  NOT 

DO  21  J= 1 ) N 
DO  21  K=1,N 
21  Pit J,<)-P(J,K> 


C  CALCULM  i£  AND  PRINT  J(R),  aI\)  t  AND 


DO  23  I = 1 , N 
23  AT  C 1*1) =ALr AT (  1,1) 

DO  28  J  =  2  » N 
DO  28  K=1 »N 
28  AT(J,<)=0.0 

CALL  PR0DIU»AT.X»1.1.N) 
Z=U {  1 , 1  ) 

CALL  LOST  ( DuL ,X»0»R,Z,N) 
1  COS  T  =  TCO  S  T  +  DOL 

L=  1 


A  N  D  P  ( .<  )  . 


N  ) 


COST  tC). 


APPENDIX  III 


I  I  1-3 


PR OGR A. -i  0PC0N2  (  COST  I  No't_D  ) 


8? 


25 


26 


29 


27 


2  4 
1 
2 
5 

4 

5 

6 
7 
6 
9 


PRINT  7»L,Z.(X(I»1)»I=1,N) 

PRINT  9  ( T COST ) 

KN=1 

DO  24  1=2) M 
DO  2  5  J  =  1 , N 
Y  (  J»  1) =X( J  > 1 ) 

CALL  PROD (X. PH  I *Y*N» 1 *N) 

DO  26  J=1,N 

DEL  T  A ( J  ,  ]  }  =  Z*DEL ( J.  I  ) 

X  t  J )  1 } =  X( J) 1 ) +  D F L T A ( J  ,  1  } 

AT  (  1 , J } =AL F AT t I » J) 

DO  29  J=2»N 
DO  29  K  =  1  ) N 
AT ( J*K ) =0.0 

CALL  PROD( U » AT , X, 1 , 1 *N ) 

Z=U (1)1) 

CALL  COST  (DOLiX»Q.R»Z.M 
T  COST  =  TCOST  +DOL 
AN= KN+ 1 

IF  ( KN- 11)24)27)27 

PRINT  7»I ,Z» (X( I ,1 ) » I=1,N) 

PRINT  9  (TCOST) 

KN  =  0 

CONTINUE 
FORMAT  ( 8F1 0 • 6 ) 

FORMAT  (15) 

FORMAT  (IF  10.6) 

FORMAT  ( / / ) 2  HM=  » I  5 ) / ) 2  H A  T  »8X»3rlu«6»//  ) 
FORMAT  <//.2HM=*I5./.3HP3I » 7X » SF10 . 6  ,  / /  ) 
FORMAT  ( // »ZHH= » I  5  » / > 1  HP » 9X >  6 F 10 . 6  > / f  ) 
FORMAT  (  //  »2riM=  » I  5  »  3X ,2HU=  »lrlG.6*/»iGX» 
FORMAT  (  3F10.6 ) 

FORMAT  ( SOX ) 6HTCOST= ) 1 F 15 .6 ) 

END 


3F10.6)// ) 
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II  1-4 


PROGRAM 


OP COM2  (CONTINUED) 


C  THIS  SUBROUTINE  CALCULATES  A ( K ) T . 

SUBROUTINE  ATRAN  ( AT  ,P tPHI .DEL.R.N) 

OOIHFNSTON  AT ( 8 * B  >  *  P ( 8 * 8 ) ) PH  I  ( ti , 8 ) »  DEL (8)8 ) »DElT ( 8. H » *Ao( 8, 8 ) 
1 AC ( 8 • B )  « AD ( B  »  8 ) 

CALL  TRANCOl.  (  DEL  »OE  L  T  *  N  ) 

CALL  P ROD ( Ad  t  T  L  T  »  P »  1 »  N  » N ) 

CALL  PROD ( AC i  *  DE L * i * i * N ) 

CALL  P RC D  (  A D  >  f i o  i  P n  I  *  i  » N  *  N  ) 

DO  1 A  I = 1  ,  N 

14  AH  1. I  ) =~AD (  1*1  )/( ACT  1 *1 J+R) 

END 

C  THIS  SUBROUTINE  CALCULATES  PSI(K). 

SUBROUT  I N E  PPS I  ( PS  T  *  PH l . DEL  *  A i  *  H ) 

DIMENSION  PS!  (  6  *  B  )  » P  Pi  I  (  a  *  d  }  *  D  F_  L  (  »  *  6  )  *  A  T  (  6  *  6  )  » A  B  (  3  *  6  ) 

TALL  PROD  ( A  A  *  D 1  L  *  A  I  » N  *  N  *  1  > 

CALL  SUM  < PSI >PHI * AB,N ) 

END 

C  THIS  SUBROUTINE  CALCULATES  PIN). 


SUBROU I  INF  PP ( P  *  PS  I . P 1 » 0  *  A 1  ,  R  *  N  ) 

8  J  I  MENS  I  ON  P  (  3  *  6  )  *  P 1  (  b  *  8  )  » PS  I  (  6  *  6  )  *0(6*8)  » AT  (  6  •  6  )  *  AA  (  6  *  ti  )  * 

I  P  S 1  T  (  8  *  8  )  ,  AC  (-8*8)  *  AO  <  6  »  8  )  ,  AE  (  8  *  8  ! 

DO  Id  I  =  1  *  N 

1 5  AA ( I  *  1  )  =A1  (  1 , I  ) 

CA  L  L  T  RAH  SO  (  P-S  I  ,  PS  I  T  ♦  N  ) 

CALL  PROD  (  A  3  *  P  3  I  T  *  r  1  *  N  *  N  *  N  ) 

CmLL  0  A  0  D  (  A  *  A  d  *  i J S  I  *  T  *  N  *  lj) 

CALL  r-rtCDt  AJftrtw  .AT  »N  »N  ,1  )  * 

DO  16  I - 1  » N 
DO  1 o  J -  i  *  N 

16  A  D (  I  » J ) =  K  * A  D (  I  * J ) 


App- NO  I  X  III 


III  -5 


P  A  0^  ic  A  ,• 


C  ONE  (  CON  I  I  "JUEO  ) 


n  n 


CALL  SUM  (AF,AC,0,R) 
CALL  SUM t  P »  A£ » AD » N ) 
END 


C  THIS  SUBROUTINE  TRANSPOSES  A  COLUMN  MATRIX  HAVING  A  MAXIMUM  OF 

c-  eight  elements. 

SUBROUTINE  TRANCCL  ( A  » 3  ,  N  ) 

DIMFNS  ION  A  (  8,  b  )  ,  fj  {  rt  *6  I 
DO  12  I =1  *  N 
12  Bl 1  ,  I ) =A(  I . 1 ) 

END 

C  THIS  SUBROUTINE  TRANSPOSES  A  SQUARE  ■MATRIX  OF  MAXlnjM  ORDER  bXB. 

SUBROUTINE  T  RAN SO  (  A  ,  rt  ,  N  > 

DIMENSION  A(8»6)»3(3»8) 

DO  II  1=1 , N 
DO  II  J  =  1 » M 
11  6  (  J  ,  I  )  =  A  <  J  »  J  ) 

END 

C  THIS  SJJBROUT I  N.£l  CALCULATES  Tr;E  COST  AT  EACH  ITERATION  AS  DETERMINED 
C  BY  THE  COSI  FUNCTION. 

subroutine  cost  <dol,y,d,r»z,n) 

D  I  i-i  lN5  ION  Y  (  6  »  6  )  ,016,3)  )Z  I  oiu  ]  ,  i  T  (  b  »  a  )  ,  A  (  o  »  c  j  ,  mm  (  g  ,  g  ) 

CALL  TRANCOL (Y.YTiN) 

CALL  PRuD (A,YT»Q»1,N»N) 

CALL  PROD  I AA,A»Y » 1 • 1 ,N ) 

DO  L  =  A  A ( 1,1) +  R #Z*l 
END 

THIS  3  U  B  R  0  u  i  I  N  z.  Ri  J  L  T  I  P  L  I  c  3  ANY  T  .vO  MA  T  R  I  C  c  o  vv  H I  C  n  /-IRE  Li  M I  T  E  D  TO 
TWO  8X6S.  Thu  ARGUMENTS  ARE  DEFINED  AS  FOLLOWS, 


APPENDIX  III 


m-6 


PROGRAM  OP COM2  (CO^TInOcD) 


n  r\  n 


A  T  Hr  PRODUCT 
5  THE  MULTIPLICAND 
C  THE  MULTIPLIER 

L  NUMBER  OF  ROWS  OF  THE  MULTIPLICAND  AND  PRODUCT 

'•I  Number  of  columns  of  the  multiplier  and  product 

N  NUMBER  Or  COLUMNS  OF  The  MULTIPLICAND  A N D  Trie  NUMdER  uF  ROWS 
OF  TnE  MULTIPLIER 

SUBROUTINE  PROD  (Atd»C»L»M»N) 

D I M E  NS  I  ON  A(6»81*s(3»d)»C(6»6) 

DO  10  I =1 »  L 
DO  10  J=1,M 
A  (  I  j  J  I  =  0  ■  0 
DO  10  K  =  1  *«4 

1U  A(  I  .  J1=A<  I  »  J)+=>(  I  .M*CCK-«J) 

END 

C  THIS  SUBROUTINE  FINDS  TrlE  SUM  OF  ANY  TWO  SQUARE  MATRICES  OF  .THE  SAME 
§5  C  D I  Fi c, R o  IONS  UP  "i  G  AH  oXB  • 


C  (A) 
C  {  H) 
C  (  C) 
C  (  D) 
{  E ) 
(  F  ) 


SUBROUTINE  SUM  IA»B»C»N) 

D I  MENS  1  OH  A ( a >  6 ) »  6 ( d *  a ) »  C ( o  *  d } 
DO  13  I =] »N 
DO  13  J  =  1  *  N 

13  A(  I  ,  J  )  =RU  ,  J  n-C  (t  .  J  > 

END 

END 
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nr  i 

(THE  SYSTEM  IMPULSE  RESPONSE  MATRIX,  FORMAT 
IN  BY  COLUMN) 

3  c  1  3 . 6  ■  READ 

r 

w 

r 

(  F) 

PHI 

(THE  SYSTEM  TRANSITION  MATRIX,  FORMAT. NF10. 
ROWS  ) 

6,  READ  IN  BY 

f 

(  G) 

0 

(WEIGHTING  MATRIX  USED  IN  THE  COST  FUNCTION, 

FORMAT  MF10.S 

C  READ  IN  3Y  ROWS) 

C  FORMAT  STATEMENT  MR  5  MILL  HAVE  TO  BE  ALTERED  DEPENDING  UPON  THE 
C  0QDER  OF  THE  SYSTEM.  I T  ■=■  FORM  ACT  SHOULD  BE  NF10.6. 


C  READ  DATA  AMD  INITIALIZE  VARIABLES. 

READ  2  (  ■  i) 

READ  2  IN) 

READ  2  (IT) 

READ  3  { R ) 

RE^D  1  (DELI  I  »  1  )  •  I  =  1  ,  N  ) 

READ  8  ( ( PH  I (  I  ,  J }  *  J~ 1 1 N )  ■  I  -  1  •  N ) 

READ  3  (  (  0  (  I  .  J  }  .  J  -  1  ,  M  )  •  I  ="1  ,  "  ) 

PRINT  2  CM) 

PPJNT  7  ( N ) 

PRINT  2  ( IT ) 

PRINT  3  ( R ) 

PRINT  ]  ( D*L( I  -  1 ) ? 1  =  1 , N) 

PRINT  S  ( (PHI ( I . J) , J  =  1 ,N)  j  I  =  1,N) 

PRINT  8  ( (Qf I , J  J , J  =  1 »NI f I =1 »M ) 

KN  =  0 

DO  2  0  I = 1  ,  N 
DO  20  J  =  1  ,N 
20  PI ( I f J ) =0( I » J) 

C  CALCULATE  AND  PRINT  A  (  K  )  T  ?  PS[(<)  .  AND  p  (  ■<  j  , 
DO  21  1=1, M-l 


APPENDIX  IV 


I  V-2 


PROGRAM  OPCON3  (CONTINUED) 


CALL.  ATRANl  A  T  *P1  ,  PHI  ,  DEL,R  ,M  ) 

CALL  P^SI  (  nSI  ,PHJ  •&€■(_»  U  »N  ) 

CALL  PD(P»PSI>°lfp.AT«9jV) 

KN=KM+1 

IT  t KM- IT ) 17,18,18 
18  PRINT  4  » I , ( AT ( 1 , J  ) » J  =  1 ,N) 

PRINT  5*1,1  (  PS  I  (  J,  K)  ,K>1  >N  )  ,  J  =  1  ,  N  ) 
PRINT  6  »  I  ,  (  !  P  (  J  ,  K  )  j  K.  =  1  ,  N  )  »  J- 1  *  R  ) 
KN-0 


17 

CO NT INU 

r 

DO  21  J 

=  1  *  N 

DO  21  K 

=  1 ,  N 

21 

PI ( J»X) 

=  P(  J  , 

K  ) 

1 

FORMAT 

{ 4  F 1 5 

.10) 

2 

FORMAT 

(IS) 

■3 

format 

!  1  c‘  1  n 

.5  ) 

4  FORH  AT  t//,2HM=-I2,/r2HAT,8X,?c‘l,''.6-//1 

5  FORT  AT  (  //  •  2IJ>'-  •  12*  /  ,  3HPS1  WX  »  6cn.  6,  //  ) 
5  FORMAT  (  //  ,  2HM=  ,  I  2  ,  /  ,  1U  J  -  NX  -  8  FID . 5  -  /  /  ) 

8  FORMAT  f  2 F 1 0*6 ) 

END 

C  THIS  SUBROUTINE  CALCULATES  At  KIT. 


SUBROUTINE  ATP  AN  (AT  *  P , °H I  •  DC L , R , N  ) 

0  R  iv  c  NS  I  ON  A  r  (  3 . 3  )  . r’  [  r  .  a  }  ,  pH  i  (  q  .  r>  i  ,  PEL  !  R  ?  8  !  ,  0  EL  T  (  8 , 8  }  ,  ' 
!  *• C  (  8  •  3  )  •  '  R  {  8  *  °  ) 

CALL  TRAM COL  ( DEL r PE LT , M ) 

CALL  PROD ( AB • DELT • P , 1,V,N) 

CALL  PROD ( AC  »  4U , PEL  > 1 >1 s  N) 

CALL  PROD  (  *  P ,  A  D' ,  PH  I  >  1 ,  N  *  N  ) 

DC  14  I = 1 , N 

14  AT { 1 , I  ) =- AD { 1 • I  ) / ( AC ( 1 . 1  I +R I 
END 
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